#include <x86intrin.h>
#include <stddef.h>

#define UNROLL 4
void dgemm_avx_unroll(size_t n, double *A, double *B, double *C)
{
	__m256d c[UNROLL];
	for (size_t i = 0; i < n; ++i) {
		for (size_t j = 0; j < n; j += 4 * UNROLL) {
			c[0] = _mm256_load_pd(C + i * n + 4 * 0 + j);
			c[1] = _mm256_load_pd(C + i * n + 4 * 1 + j);
			c[2] = _mm256_load_pd(C + i * n + 4 * 2 + j);
			c[3] = _mm256_load_pd(C + i * n + 4 * 3 + j);

			for (size_t k = 0; k < n; k++) {
				__m256d a = _mm256_broadcast_sd(A + i * n + k);
				c[0] = _mm256_add_pd(c[0],
						_mm256_mul_pd(_mm256_load_pd(B + k * n + 4 * 0 + j),
							a));
				c[1] = _mm256_add_pd(c[1],
						_mm256_mul_pd(_mm256_load_pd(B + k * n + 4 * 1 + j),
							a));
				c[2] = _mm256_add_pd(c[2],
						_mm256_mul_pd(_mm256_load_pd(B + k * n + 4 * 2 + j),
							a));
				c[3] = _mm256_add_pd(c[3],
						_mm256_mul_pd(_mm256_load_pd(B + k * n + 4 * 3 + j),
							a));
			}

			_mm256_store_pd(C + i * n + 4 * 0 + j, c[0]);
			_mm256_store_pd(C + i * n + 4 * 1 + j, c[1]);
			_mm256_store_pd(C + i * n + 4 * 2 + j, c[2]);
			_mm256_store_pd(C + i * n + 4 * 3 + j, c[3]);
		}
	}
}
